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Abstract: We propose a new Monte Carlo technique in which the de- 
generacy of energy states is obtained with a Markovian process analogous 
to that of Metropolis used currently in canonical simulations. The obtained 
histograms are much broader than those of the canonical histogram technique 
studied by Ferrenberg and Swendsen. Thus we can reliably reconstruct ther- 
modynamic functions over a much larger temperature scale also away from 
the critical point. We show for the two-dimensional Ising model how our new 
method reproduces exact results more accurately and using less computer 
time than the conventional histogram method. We also show data in three 
dimensions for the Ising ferromagnet and the Edwards Anderson spin glass. 

In the simulation of thermal systems one typically needs to calculate a variable < Q >t 
as a function of temperature T. Usually one would have to perform independent Monte 
Carlo simulations at different values of temperature. An appealing strategy to avoid such 
multiple simulations is the canonical histogram method [1,2]: The histogram Pt(E) of the 
energies at one given temperature To is measured and then the distribution at a different 
temperature T is obtained by reweighting, i.e. by multiplying with exp(E/T — E/Tq) and 
normalizing. In order to obtain the average < Q >t, one needs to accumulate also the 
histogram of Q as a function of energy E. The thermal average at temperature T is then 
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where < Q{E) > means the average value of Q obtained at fixed energy E. For simplicity 
we set ks = 1 in this article. 

The histogram method has been carefully checked on various models [3]. Its main 
disadvantage is that the canonical distribution Pt(E) of the energy is rather narrowly 
peaked around the average value < E >t q (and the more so the larger the system) so 
that when T is not close to Tq the reweighting factor is very small there and very large 
near < E >y. One also has strong statistical fluctuations stemming from the tails of the 
distributions. So, in order to get reliable results the tails must be sampled very well by 
making good statistics and the temperatures considered should not be far from T . 

We want to introduce here a completely different technique based on the calculation of 
the degeneracy g(E) of energy states from histograms of adequate macroscopic quantities 
defined below [4]. Figure 1 compares the classical canonical histogram for an Ising model on 



square lattices to the one obtained with our method. As mentioned before the width of the 
classical histogram method decreases with system size. Instead the histogram introduced 
in the following covers the entire energy range. 

The first step of our method [4] is to perform a Markovian walk along the energy axis 
which samples all energies with the same weight. For that we define two classes of moves 
in phase space: 

class 1: E ► E - AE 

class 2: E ► E + AE , 

where AE > 0. To correctly sample g(E) we propose the following dynamics: We randomly 
choose a move. If this move belongs to class 1, it is accepted. If, however, it belongs to 
class 2, it is accepted only with probability N dn /N up , where N dn and N up are the total 
numbers of possible moves of classes 1 and 2, respectively, measured at the current state. 

These quantities iVdn and N up concern the whole lattice, considering all possible move- 
ments one can perform at the current state in order to get a (would-be) next state along 
the Markovian chain. In the particular case of a single spin flip dynamics, N^ n and N up are 
obtained by checking on each site if a spin flip would increase or decrease the energy and 
counting how many times each case occured. Then one site is chosen randomly and the 
above rule is applied only to it. Nevertheless, the method is not restricted to single spin 
flips, and any other updating protocol can be adopted, provided all possible movements 
are classified as belonging to classes 1 or 2, in order to count their current numbers N dn 
and iV U p, before actually performing one of them. Following this rule the probabilities to 
increase or to decrease the energy are equal. 

Like for a random walk, the region already visited along the energy axis increases 
proportionally to AEy/i, where t is the number of performed moves, i.e. the length of the 
Markovian sequence of states. This allows one to obtain any predefined distribution width, 
simply by taking a large enough computer time. On the other hand, in canonical dynamics 
based on the Boltzmann weights one always get distributions with exponential decaying 
tails, and thus with time-independent widths. Figure 1 shows the number of visits as a 
function of the energy. 

The second step of our method is the calculation of the degeneracy of energy states 
g(E). In the above Markovian process, the probability for the energy to jump from E to 
E + AE is the same as that of jumping back from E + AE to E: 

< N up (E) > g(E) = < N dn (E + AE) > g(E + AE) (2) 
where the averages are taken at fixed energy. Equation (2) can be rewritten as 

_ d\ng(E) _ 1 < N up (E) > 
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In this way one obtains g{E), from the averages < N up (E) > and < Nd n (E) > accumulated 
during the random walk. 

Summarizing, the method consists in performing the Markovian process defined above, 
and accumulating four histograms along the E axis: for the number of visits; for the 
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quantity Q one is interested in; for the average number A^n of moves of class 1; and for 
iV U p corresponding to class 2. From these numbers, g(E) is determined by equation (3), 
and the thermal average < Q >t obtained through equation (1) using 

P T (E) = ^g(E)eM-E/T) (4) 

with Zt = ^ E g{E)ex.\)(—E/T), if one wants to work in the canonical ensemble. 

For the one-dimensional Ising model we have been able to prove that our model gives 
the exact solution for g(E) [4]. Also in this same reference one can find a generalization 
of equation (3) for the case where different values of AE are possible: it is enough to 
transfer the factor 1/AE as an exponent in both the numerator and denominator, inside 
the average brackets. 

In order to test our new technique numerically we performed simulations for three 
examples: the two- and three-dimensional Ising ferromagnet and the 3D Ising spin glass. 
Starting with a random configuration we choose in the ferromagnetic cases an initial state 
that was thermalized at the critical temperature for 10 Metropolis lattice sweeps [5]. For 
the spin glass we did instead 100 sweeps at zero temperature. We checked, however, 
that the results did not depend on the initial state. Our systems had periodic boundary 
conditions in all directions. For efficiency we used the multilattice approach [6] storing 32 
lattices in an array of 32-bit integers. All 32 samples are processed in parallel by using 
multi-spin coding [7]. 

The degeneracy g(E) of states obtained for the two-dimensional Ising model on a 
32 x 32 square lattice agrees so well with the recently derived exact degeneracy [8] that 
on a plot they are indistinguishable from each other. Figure 2 shows the averaged energy 
and specific heat, obtained by the present method. The data are extremely close to the 
exact curves [9]. This is true even for the specific heat at its maximum where statistical 
fluctuations are largest as can be seen from the inset of figure 2. 

In figure 2 we also see data obtained using the classical canonical histogram method 
[1,2]. The data noticeably deviate around the maximum of the specific heat which is 
probably due to the fact that the width of the canonical histogram is not larger than 
the distance between this maximum and T c , the temperature at which the histogram was 
calculated. Although the computational effort invested in CPU time was about 50% larger 
for our method than for the canonical histogram method (using in both cases the same 
multi-spin implementation) the result obtained with our method is at least ten times more 
precise. Moreover, for the canonical histogram technique one must simulate independently 
at several temperatures to get the data in the temperature range considered. This is a 
convincing evidence for the efficiency of our method. 

For the two-dimensional Ising ferromagnet the magnetization [10] and susceptibility 
were also analyzed and compared with canonical Monte Carlo simulations [11]. The agree- 
ment is again excellent. We did the same studies also for the three-dimensional Ising 
model. Again energy, specific heat, magnetization and susceptibility were calculated and 
found in very good agreement with known results investing a quite modest computational 
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effort. In figure 3 we show as an example magnetization and susceptibility for 3D systems 
differing in size. For large lattices like the 50 x 50 x 50 presented in figure 3, the random 
walk dynamics may be very slow, if one starts always from only one energy. In these cases, 
it is better to restart the dynamics many times during the simulation, each time choosing 
a new value for the starting energy. This trick was adopted also for the 256 x 256 lattice 
in figure 1. 

As a more difficult test we chose also to simulate the three-dimensional Ising spin 
glass. Besides energy and specific heat we also calculated the overlap order parameter 
defined through 



where {s\ } is a ground state previously determined by the approximation introduced in 
[12]. In figure 4 we show that it goes to zero close to the predicted temperature [13]. Again 
the method turns out to be very efficient. 

Two last points deserve comments. First, there are other methods [14] which can be 
thought as similar to our first step, where we defined a non-biased random walk dynamics. 
All of them introduce non-canonical ensembles directly weighted by the unknown entropy, 
or g(E). Karliner et al [14], for instance, use the g(E) scaled up from exact calculations 
on small lattices in order to define the statistical weight of their dynamics, replacing the 
Boltzmann factors. In all these earlier works, the entropy or g(E) is measured directly from 
the accumulated number of visits to each energy channel, within the particular dynamics 
corresponding to each case. Our approach is distinct from those others: we do not use the 
number of visits, on the contrary, our results are extracted from accumulated histograms 
of quantities measured within each visited state, namely N^n and A up , concerning the 
whole lattice. The only role played by the number of visits to each energy channel is that 
it must be large enough to provide a good statistics. Being based on statistical averages 
of macroscopic quantities, our method also has a better numerical precision. Indeed, 
accumulating the macroscopic quantities Ad n and A up (both of the order of the number 
of sites), instead of merely incrementing by unity a visit counter after each movement as 
in the traditional histogram method, we succeed in getting better numerical accuracy. 

Second, the left hand side of equation (3) coincides with the formal definition of 
(inverse) temperature as a function of the energy E. Thus, the averages on the right hand 
side of the same equation correspond to the canonical ensemble equilibrium. Although in 
our second step we have derived equation (3) from the particular dynamics defined in the 
first step (the non-biased random walk), it holds for any other dynamics respecting the 
canonical equilibrium for each (inverse) temperature (3(E). In this sense, the importance 
of our method relies on the second step, the determination of g(E) from the countings of 
potential movements Ad n and A up , measured within any canonical equilibrium dynamics. 

We have presented a new histogram Monte Carlo method which as compared to the 
traditional one based on temperature [1,2] is based on the histograms of the energy. These 
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histograms have the advantage of having much broader tails allowing to extrapolate to 
a much larger range of temperatures with a rather small number of samples. We have 
tested our method on the two- and three-dimensional Ising model and the 3D Ising spin 
glass and succeeded in reproducing thermodynamic quantities with higher accuracy over a 
larger temperature scale and with less computational effort than the canonical histogram 
method. 

We thank L. de Arcangelis, D. Stauffer, S. Moss de Oliveira, C. Moukarzel and Yi- 
Cheng Zhang for helpful discussions. 
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Figure Captions 



Figure 1 Number of visits as a function of energy obtained from the traditional histogram 
method [1,2] (fixing the temperature at the critical value), and from the present 
method, for the Ising ferromagnet on 32 x 32 and 256 x 256 square lattices. The 
energy is scaled such that the critical value corresponds to 0.146. Energies above 
~ 0.4 are not sampled for the 256 x 256 lattice, because they do not contribute to 
thermal averages at the interesting range of temperatures (0 < T < 2T C ). 

Figure 2 Averaged energy and specific heat obtained from the present method for the Ising 
ferromagnet on a 32 x 32 square lattice and the exact ones [8] (both full lines). Except 
for a small region around the critical point (see inset for maximum of specific heat) 
the two curves are indistinguishable. Results of similar quality were also obtained 
for larger lattices (we tested L = 64, 128 and 256). The open circles represent data 
obtained using the canonical histogram method [1,2] for which Pt(E) was obtained 
at T c . Our method needed just 40 minutes on a Pentium (66MHz) for the 6 x 10 4 
lattice sweeps performed. 

Figure 3 Magnetic susceptibility obtained from the present method for the 3D Ising ferromagnet 
on a 50 x 50 x 50 (upper curves, showing also the magnetization), a 20 x 20 x 20 (middle 
curve) and a 10 x 10 x 10 (lower curve) cubic lattice. 

Figure 4 Overlap with one specific previously calculated low temperature state and non-linear 
susceptibility of this order parameter for the 3D Ising spin glass simulated on a 8 x 8 x 8 
cubic lattice making 53000 lattice sweeps for 128 different configurations of couplings. 
The total amount of CPU time spent to get these data was 15 hrs on a Pentium 
(66MHz). 
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